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Abstract 

We present an innovative numerical approach for 
setting highly accurate nonlocal boundary condi- 
tions at the external computational boundaries when 
calculating three-dimensional compressible viscous 
flows over finite bodies. The approach is based 
on application of the difference potentials method 
by V. S. Ryaben’kii and extends our previous tech- 
nique developed for the two-dimensional case. The 
new boundary conditions methodology has been suc- 
cessfully combined with the NASA-developed code 
TLIS3D and used for the analysis of wing-shaped con- 
figurations in subsonic and transonic flow regimes. 
As demonstrated by the computational experimentsr 
the improved external boundary conditions allow 
one to greatly reduce the size of the computational 
domain while still maintaining high accuracy of the 
numerical solution. MoreoverlThey may provide for 
a noticeable speedup of convergence of the multigrid 
iterations. 

1 Introduction 

Preliminaries. External flows over finite bodies 
or configurations of bodies represent a wide class of 
important practical applications in fluid dynamics. 
To treat this type of problems numericallyEone typ- 
ically truncates the original infinite domain. The 
resulting truncated problem is obviously subdehnite 
unless supplemented by the proper closing procedure 
at the external computational boundary. The lat- 
ter procedure is called the artificial boundary condi- 
tions (ABC’s). In the ideal caselThe ABC’s would 
be specified so that the solution on the truncated 
domain coincides with the corresponding fragment 
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of the original infinite-domain solution. The issue of 
ABC’s is significant in CFD and in other areas of sci- 
entific computing; theoretical estimates and compu- 
tational experiments by different authors show that 
the proper treatment of external boundaries has a 
profound impact on the overall performance of nu- 
merical algorithms and interpretation of the results. 

Different ABC’s methodologies have been studied 
extensively over the recent two decades. Howeverr 
the construction of the ideal (i.eTexact) ABC’s that 
would provide no error associated with the domain 
truncation and at the same time be computationally 
inexpensiver easy to implementT and geometrically 
universair still remains a fairly remote possibility. 
Among the variety of approaches proposed to date 
only a few can be regarded as the commonly used 
tools in CFD. As a rulerthese approaches are based 
on the essential model simplifications (e.gT locally 
one-dimensional treatment near the external bound- 
ary) andr thereforer often lack accuracy in compu- 
tations. Thisrin turnFnecessitates choosing exces- 
sively large computational domains. On the other 
handrthese simple methods usually provide for lo- 
cal ABC ’sr andr thereforer for cheapT geometrically 
universairand algorithmically simple numerical pro- 
ceduresrwhich are attractive for practical use. 

There arer of courser methods of another kindr 
which typically provide for highly accurate and ro- 
bust numerical algorithms. These methodsr how- 
everrare not used routinely because the correspond- 
ing ABC’s in most cases appear global. As a conse- 
quencerthese ABC’s may be relatively cumbersome 
and expensive; moreoverrthey can be derived easily 
only for the boundaries of regular shape. 

The nonlocality is inherent for exact ABC’s; on 
the level of PDE’s such boundary conditions typi- 
cally involve pseudo differential operators. As con- 
cerns the approximate approachesr the basic trend 
is the following: higher accuracy for the boundary 
procedure requires more of the nonlocal nature of 
the ABC’s to be somehow taken into account. A 
survey of methods for setting the ABC’s in different 
areas of scientific computing can be found in our re- 
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cent work 3 Tas well as in the comprehensive reviews 
by Givoli 2 - 3 . 

In this paperrwe concentrate on constructing the 
ABC’s for the three-dimensional problems of compu- 
tational aerodynamics. Werhoweverr mention that 
this area constitutes a fraction of the possible range 
of applications for the different ABC’s techniques. 
Besides the hydro- and aerodynamic problems (ex- 
ternal flowsr duct flowsr boundary layersr free sur- 
facesr etc.)T the entire range includes the flows in 
porous mediarhltrationrMHD flowsr plasma (e.g.T 
solar wind)rthe problems of solid mechanics (in par- 
ticularrelasticity and aeroelasticity)rand the prob- 
lems of wave propagation (electromagneticlacousticr 
seismic)rjust to name a few. 

Main Objective and Method of Approach. 

For external flow computations! 1 our goal is to de- 
rive and implement the ABC’s that would com- 
bine the advantages relevant to both global and 
local methodologies. Our approach to construct- 
ing the ABC’s is based on usage of the finite- 
difference analogues to Calderon’s generalized po- 
tentials and boundary projections (see the original 
work by Calderon ^ and also work by Seeley 3 ). 
Ryaben’kii (see Refs. 6T7F8) had modified the orig- 
inal Calderon’s construction and proposed a numer- 
ical technique for the effective calculation of poten- 
tials and projections; this technique is known as the 
difference potentials method (DPM). 

In Ref. 9 we describe the foundations of the 
DPM-based approach to setting the ABC’s for com- 
putation of two-dimensional external viscous flows 
(Navier-Stokes equations). In Ref. 10 we imple- 
ment this approach along with the multigrid Navier- 

11 19 i o 

Stokes algorithm by Swanson and Turkel ’ ’ 

and present some numerical results for subsonic and 
transonic laminar flows over single-element airfoils. 
In Ref. 14 we show the results of subsequent numer- 
ical experiments and propose an approximate treat- 
ment of turbulence in the far Held. Our work Ref. 15 
delineates the algorithm for solving one- dimensional 
systems of ordinary difference equations that arise 
when calculating difference potentials and the DPM- 
based ABC’s. In Ref. 16 we extend the area of appli- 
cations for the DPM-based ABC’s by analyzing two- 
dimensional flows that oscillate in time; we also pro- 
vide some solvability results for the linearized thin- 
layer equations used for constructing the ABC’s. In 
Ref. 17 we present a general survey of the DPM- 
based methodology as applied to solving external 
problems in CFDrincluding parallel implementation 
of the algorithmFcombined implementation of non- 
local ABC’s with multigridrand entry- wise interpo- 
lation of the matrices of boundary operators with 


respect to the Mach number and the angle of at- 
tack. Additionallyrin Ref. 17 one can find some new 
theoretical results on the computation of general- 
ized potentialsrthe construction of ABC’s based on 
the direct implementation of boundary projections 
(thin-layer equationsjr and some numerical results 
for various airfoil flows: laminar and turbulentrtran- 
sonic and subsonicr including very low Mach num- 
bers (incompressible limit). Finallyrin Refs. 18T19 
we outline basic elements of the DPM-based ABC’s 
for the case of three-dimensional steady-state exter- 
nal viscous flows; this case is undoubtedly the one 
most demanded by the current computational prac- 
tice. Specificallyr in Refs. 18r 19 we address the 
flows around wing-shaped configurations and show 
some preliminary numerical results for the subsonic 
regime. 

Belowr we report on the subsequent develop- 
ment of the DPM-based approach for the three- 
dimensional external flows. We present the results of 
computations for different flow regimes and compare 
the performance of the DPM-based ABC’s with the 
performance of the standard boundary conditions. 
The assessment is conducted with respect to both 
the quality of the resulting solution (accuracy as it 
depends on the size of the computational domain) 
and the properties of the numerical procedure as 
influenced by the type of ABC’s (efficiencyFrobust- 
nessFand convergence of the pseudo-time iterations). 
The results demonstrate the superiority of the DPM- 
based ABC’s over the standard existing methods. 

2 External Flow Problem 

Formulation. We consider a steady-state flow of 
the viscous compressible perfect gas past a three- 
dimensional wing. The flow is uniform and sub- 
sonic at infinity; it is symmetric with respect to the 
Cartesian plane z = OTwhichrin particularrimplies 
that the free-stream velocity vector is parallel to this 
plane. 

The flow equations are integrated on the grid gen- 
erated around the wing. This grid actually defines 
the finite computational domain; the ABC’s that 
would close the truncated problem should be set at 
the external coordinate surface of the grid. Let us 
designate this surface T ; for one-block curvilinear 
boundary-fitted grid around the ONERA M6 wing 
the schematic geometric setup is shown in Figure 1. 

The outermost coordinate surface of the grid is 
designated Ti (see Figure 1); it represents the ghost 
nodes (or ghost cells for the finite-volume formula- 
tion). ClearlyFwhen the stencil of the scheme used 
inside the computational domain is applied to any 
node from TT it generally requires some ghost cell 
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Figure 1. Schematic geometric setup for the three-dimensional case. The wing on the left is enlarged. 


data. Unless the required data are providedr the 
finite-difference system solved inside the computa- 
tional domain appears subdefinite (i.e.Tit has less 
equations than unknowns). From the viewpoint of 
what the solution technique isF one can say that 
when some iterative solver is employed to integrate 
the flow equations inside the computational domainT 
the values of the solution at the ghost cells should 
be prescribed at each iteration in order to be able to 
advance the next “time” step. 


the possibility of linearization in the far field re- 
quiresFgenerally speakingFsome additional analysis. 
For the model full-potential formulationFwe present 
here a simple argument that the compressible three- 
dimensional far field remains linear for transonic flow 
regimes as well. 

Consider the Karman^Guderley equation (see ^°) 

d 2 <f> d 2 <j> d 2 <j> n + ld<f>d 2 <f> . 

dx 2 dy 2 dz 2 K dx dx 2 


ThereforeFin the practical framework the closure 
of the discretized truncated problem means the spec- 
ification of the solution values at the ghost cells. 
This will be done by means of the DPM-based ABC’s 
so that the boundary data used for the closure ad- 
mit an exterior complement that solves the problem 
outside the computational domain (see below). 

FirstF we assume that the flow perturbations 
against the constant fn > stream background are 
small in the far field and consider the accordingly 
simplified problem outside the computational do- 
main (i.e.routside F). ClearlyFfor the small to mod- 
erate free-stream Mach numbers Mo (purely sub- 
sonic flows) we can retain only the first-order terms 
with respect to perturbations in the governing equa- 
tions. Thereforer the aforementioned simplification 
of the original problem will actually be the far- • 
field linearization. Howeverrin the transonic limit 


for the perturbation <f> the full potential $ of 
the compressible gas flow around a thin three- 
dimensional wing. Here 


« 0 


«0 



V = s 1/3 y, y = S 1/3 y , 
y = \fKy, z = 


1 _ Mo 2 
V £2/3 ’ 

\fKz, 


( 2 ) 


8 is the wing thickness (£ — +0 along with Mo — *■ 
1 in the transonic limit )TK is the parameter of tran- 
sonic similarity (the true linear theory corresponds 
to big values of K )F // : is the flow velocity at infin- 
ityF k is the ratio of specific heatsF and the positive 
x direction coincides with the free stream. 

on 

The common practice ■ of developing the asymp- 
totic expansions for solutions of the equations that 
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involve transonic nonlinearities consists of substi- 
tuting the leading linear term(s) into the nonlinear 
parts of the equation (right-hand side of ( 1 )) and ob- 
taining the corresponding corrections by solving the 
resulting non-homogeneous problem. In the specific 
case under studyrthe linear far-held expansion starts 
with the horseshoe vortex 

< f >1 = ~TV ~2 + r = \/x 2 + y 2 + z 2 . (3a) 

Substituting the expression (3a) into the right-hand 
side of (1) and solving the resulting Poisson equa- 
tion (this involves the Fourier transform in spherical 
functions) one obtains the nonlinear correction 

<f>iNL ~ r~ 3 (3b) 


that decays at infinity two orders of magnitude faster 
than the source term (3a). Analogouslyrfor the gen- 
eral doublet term <f > 2 ~ f~ 2 the corresponding non- 
linear correction can be shown to be <f> 2 NL ~ r -5 r 
which is the three orders of magnitude difference. 
We therefore conclude that the transonic nonlin- 
ear corrections can be neglected when analyzing the 
compressible far Held in three space dimensions. 

Of courser having established the fact of the far- 
field linearityFwe cannot say in advance whether or 
not the linearization outside T is possible for every 
specific configuration of the domains. Clearlyr for 
the very large computational domain one can lin- 
earize the flow outside TT and as we approach the 
wing (the source of perturbations)r the validity of 
linearization can be verified a posteriori (seer e.g.T 
our previous work 14. If) 

The linearized dimensionless thin-layer equations 
can be written as 

dp du dv dw ' 

dx dx dy dz 

du dp 1 f d 2 u d 2 u\ 

dx dx Re \dy 2 dz 2 ) 

dv dp 1 / 4 d 2 v 

dx + dy Re \3 dy 2 


a v 
dz 2 


dydz 


dw 

dx 


dp 

dz 


1 

Re 


3 dz 2 


dy 2 


dp 

1 dp 

K 

dx 

M 2 dx 

Re Pr 

d 2 p'' 

\ 1 

(d^p 

dz 2 , 

1 kM 2 

\dy 2 


dydz 
d 2 p 
dy 2 
d 2 p 
dz 2 


= 0 


(4a) 


= 0 


= 0 


where pri/TiTr/Tand p are the perturbations of den- 
sityr Cartesian velocity componentsr and pressurer 
respectivelyrffe is the Reynolds numberrand Pr is 
the Prandtl number. System (4a) is supplemented 
by the homogeneous boundary condition at infinity: 

u =(p,u,v,w,p) — *■ (0,0, 0,0,0) (4b) 

as (x 2 + y 2 + z 2 ) — >■ 00 , 

which corresponds to the free stream limit of the 
solution. Returning to the question of closing the 
discretized truncated systemr we clarify that the 
boundary data provided by the DPM-based ABC’s 
will admit an exterior complement that would solve 
the discrete counterpart of (4a) and meet boundary 
condition (4b) in a certain asymptotic sense. 

DPM-based ABC’s — Main Idea. We dis- 
cretize system (4a) on the auxiliary Cartesian grid 
with the second order of accuracy; we use first-order 
differences in x and second-order differences in y 
and z (see Refs. 16T 17 for the details in two di- 
mensions). The DPM will provide us with the com- 
plete boundary classification (in terms of the appro- 
priate traces) of all those and only those exterior 
grid vector-functions that solve the discrete coun- 
terpart of (4a) outside the computational domain 
and satisfy boundary condition (4b) in some approx- 
imate sense. The foregoing boundary classification 
will be obtained as an image of a special projection 
operatorPwhich can be considered as a discrete ana- 
logue of Calderon’s pseudodifferential boundary pro- 
jection As we solve the Navier-Stokes equations 
inside the computational domain iterativelyr every 
time we need to update the ghost cells we take a 
certain sufficient set of data from inside T (see be- 
lowjTproject it onto the right manifoldri.eTonto the 
subspace in the entire space of boundary data that 
admit the correct exterior complementT and obtain 
the ghost cell values by calculating the trace of this 
complement on Ti . 

DPM-based ABC’s — Specific Implementa- 
tion. Let us split the nodes of the auxiliary Carte- 
sian grid into two distinct groups: those that are 
inside T and those that are outside. Applying the 
stencil of the scheme for (4a) to each node of both 
groupsrwe consider the intersection of the grid sets 
swept by the stencil. This intersection is called the 
grid boundary 7 ; it is a multi-layered fringe of nodes 
of the auxiliary Cartesian grid located near T. Fig- 
ure 2 shows an example of the grid boundary 7 (sev- 
eral cross-sections in different directions). 
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Figure 2. Grid sets for the three-dimensional case. 


For any function u on the Cartesian grid we de- 
fine its trace Tr 7 u on 7 as merely a contraction. 
For any grid function u 7 specified on 7 we introduce 
the generalized potential P u 7 with the density u 7 ; 
the generalized potential is defined on the auxiliary 
Cartesian grid on 7 and outside it. The generalized 
potential is obtained as a solution of the special aux- 
iliary problem (AP); solution of the AP replaces and 
extends the operation of convolution with the funda- 
mental solution in the classical potential theory. I lie 
AP is driven by the right-hand side that depends on 
u 7 rt.he formal construction of this right-hand side 
is the same in two- and three-dimensional cases and 
we refer the reader to our previous work ^ for 

details. Boundary conditions for the AP should ap- 
proximate boundary condition (4b) and at the same 
time ensure the finite-domain formulation for the 
AP. ThereforeFwe formulate the AP on a sufficiently 
large parallelepiped aligned with the Cartesian di- 
rections; the parallelepiped fully contains Ti. We 
specify the periodicity boundary conditions in the y 
and z directions and also assume symmetry at z = 0 . 
After the Fourier transform with respect to y and z 
the discrete counterpart to (4a) can be written as a 
family of one-dimensional difference equations: 


Ak u m,k + Bk u m-l,k — f m -l/2,ki (5) 

m = 1, . . . , M , k = (k y ,k z ), 
ky 0, . . . , Jy, k - 0 , 

where Ak and Bk are the 5 x 3 matrices and M + IT 
2J y + lTand J~ + 1 are the numbers of grid nodes 
in the xTyT and z directionsr respectively (symme- 
try is taken into accountT as well as the fact that u 
and f are real-valued). Boundary conditions in the 
x (stream wise) direction are specified separately for 


each pair of wavenumbers k: 
s “( k ) u o,k = °, ( 6 a) 

S+(k)u M k = 0 , ( 6 b) 

where 

S-(k)= n (Qk-Mk)I), (7a) 

|;i s ( k ) | > 1 

S+(k)= n (Qk-Mk)I), (7b) 

|;i s ( k ) | < 1 


Qk = A k lB k r and p s (k) are the eigenvalues of 
Qk The senri-analytic boundary conditions ( 6 a) 
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and (6b) (the eigenvalues for (7) are calculated nu- 
merically) explicitly prohibit growing modes of the 
solution in the left and right directionsFrespectively. 
The periods in y and z should be chosen sufficiently 
large to ensure that the periodic solution considered 
near T and Id is sufficiently close to the theoretical 
non-periodic solution; the latter can be thought of 
as a limit when the periods approach infinity. The 
approximation of a non-periodic solution by the pe- 
riodic one on a finite fixed interval as the period(s) 
increase(s) is discussed in our work 9 - 16 ' 17 . In 
Ref. 17Twe also discuss the possibility to replace the 
Fourier transforms by the non-unitary transforms. 
The latter appear when the grid in y or z is stretched 
(which provides for a drastic cost reduction) and 
the corresponding eigenfunctions consequently form 
a skew basis. 

The foregoing AP allows us to calculate the gener- 
alized difference potential P u 7 for any grid density 
u 7 specified on j. The composition of the operators 
Tr 7 and PTP 7 = Tr 7 Pris a projectionTP^ = P 7 T 
and it is a discrete counterpart of Calderon’s bound- 
ary projection for system (4a). The image of this 
projectionrimP 7 rcontains all those and only those 
u 7 ’s that are the traces of some exterior difference 
solution to (4a) that satisfies the boundary condi- 
tions of the AP (periodicity in y and z and boundary 
conditions (6) in x). The latter boundary conditionsr 
in turnrapproximate (4b). 

Ffaving constructed the procedure for calculating 
the potentials and projections for the discrete ver- 
sion of (4a)Twe can now close the system inside the 
computational domainri.e.robtain the ABC’s. FirstT 
we take u and du/dn on TTn is the normair(these 
data are available from inside the computational do- 
main) andrusing interpolation Rr along T and the 
first two terms of the Taylor expansion (denoted 7T 7 )r 
obtain u 7 : 


I- — 7r 7 Rr 


<9u 

dn 


( 8 ) 


ThenFwe need to calculate the potential P v 7 for the 
density v 7 = P 7 u 7 and interpolate it to the nodes 
Ti: 


u 


Ti 


RrqPv 7 = RrqPu 7 . 


(9) 


Finallyrthe ABC’s are obtained in the operator form 




= T 


<9u 

dn 


(10) 


where T is composed of the operations (8) and (9). 
Boundary condition (10) is applied every time we 


need to update the ghost cell values in the course of 
the iteration process. The implementation of ABC’s 
(10) can either be direct or involve preliminary cal- 
culation of the matrix T. In the latter caserthe run- 
time implementation of the ABC’s (10) is reduced 
to a matrix-vector multiplication. 

3 Numerical Experiments 

Two-Dimensional Case. For the reason of com- 
pletenessFwe first briefly comment here on the two- 
dimensional results from our previous work (see 
Refs. 10ri4ri7). In that workrwe calculated sub- 
sonic and transonic viscous flows past single-element 
airfoils (NACA0012 and RAE2822). 

The computational domain is formed by the C- 
type curvilinear grid generated around the airfoil. 
On this gridr the Navier- Stokes equations are in- 
tegrated using the code FLOMG by Swanson and 
Turkel 77, 72 ’ 73 . The standard treatment of the ex- 
ternal boundary in the code FLOMG is based on the lo- 
cally one-dimensional characteristics analysisFwhich 
may or may not be supplemented by the point- 
vortex (p.-v.) correction 27 . 

Basic conclusions that could be drawn from our 
two-dimensional numerical experience are the fol- 
lowing. The DPM-based ABC’s are geometrically 
universairalgorithmically simple and easy to imple- 
ment along with the existing solver. For the large 
computational domains (30-50 chords of the air- 
foil)rthe performance of the standard methods and 
the DPM-based ABC’s is roughly the same. Asr 
howeverrthe artificial boundary approaches the air- 
foil the discrepancy between the corresponding so- 
lutions increases. The lift and drag coefficients ob- 
tained on the basis of boundary conditions (10) de- 
viate from their asymptotic (50 chords) values much 
slighter (within fractions of one percent) than the 
coefficients obtained using local ABC’s do. In other 
wordsrthe nonlocal DPM-based ABC’s allow one to 
use much smaller computational domains (as small 
as 2-3 chords) than the standard boundary condi- 
tions do and to still maintain high accuracy of com- 
putations. Moreoverr if we compare three models: 
DPM-basedrpoint-vortexFand standard locairthen 
it turns out that the DPM-based ABC’s display the 
best performance for small computational domainsr 
the performance of the local characteristic boundary 
conditions for the small domains is very poorT and 
the point- vortex boundary conditions perform much 
better for the lift than they do for the drag coef- 
ficient. This behavior seems reasonable since the 
point-vortex model is a lift-based treatment. 

We also note that for the certain variants of 
computation the DPM-based ABC’s may noticeably 
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Table 1. ONERA M6: M 0 = 0.5; Re 0 = 11.7 • 10 6 ; a = 3.06°. 


Domain “radius” 

1.25 root chords 

2 root chords 

10 root chords 

Grid 

197 x 49 x 33 

Type of ABC’s 

standard 

DPM 

standard 

DPM 

standard 

DPM 

Full iiftrc L 

0.2218 

0.2065 

0.2185 

0.2065 

0.2081 

0.2072 

Relative error 

6.58% 

0.34% 

5.0% 

0.34% 

0% 

0% 

Full dragrCo x 100 

0.817 

0.791 

0.793 

0.791 

0.787 

0.788 

Relative error 

3.8% 

0.38% 

0.76% 

0.38% 

0% 

0% 


Table 2. ONERA M6: M 0 = 0.84; Re 0 = 11.7 • 10 6 ; a = 3.06°. 


Domain “radius” 

3 root chords 

10 root chords 

Grid 

197 x 49 x 33 

209 x 57 x 33 

Type of ABC’s 

standard 

DPM 

standard 

DPM 

Full iiftrc L 

0.298T0.004 

0.2798 

0.2805 

0.2786 

Relative error 

6.24%±1.43% 

0.43% 

0% 

0% 

Full dragrCo x 10 

0.168T0.008 

0.1537 

0.1542 

0.1531 

Relative error 

8.95%±5.19% 

0.39% 

0% 

0% 


speed up (by up to a factor of three) the conver- 
gence of the multigrid iterationsEsee Refs. 9riOIT4. 
The discussion on combined implementation of the 
DPM-based ABC’s with multigrid is contained in 
Ref. 17. 

Three-Dimensional Case. Herer we demon- 
strate some new results on computation of the ex- 
ternal flows around the ONERA M6 wing. 

We use the code TLIS3D by Vatsalet al. 22 to inte- 
grate the thin-layer equations on the curvilinear grid 
(see Figures 1 and 2) generated around the wing. 
This code is based on the central-difference finite- 
volume discretization in space with the first- and 
third-order artificial dissipation. Pseudo-time iter- 
ations are used for obtaining the steady-state solu- 
tion; the integration in time is done by the five-stage 
Runge-Kutta algorithm (with the Courant num- 
ber calculated locally) supplemented by the resid- 
ual smoothing. For the purpose of accelerating the 
convergencer the multigrid methodology is imple- 
mented; in our computations we used three subse- 
quent grid levels with V cycles; the full multigrid 
methodology (FMG) could be employed as well. In 
additionr we use the preconditioning technique of 
Ref. 23 to improve the convergence to steady state. 
We implement the DPM-based ABC’s on the finest 
level of multigrid on the final FMG stage; the bound- 
ary data for coarser levels are provided by the coars- 
ening procedure. MoreoverFeven on the finest level 


we implement the DPM-based ABC’s only on the 
first and the last Runge-Kutta stagesr which seems 
to make very little difference compared to the im- 
plementation on all five stages; the boundary data 
for the three intermediate stages are provided from 
the DPM-based ABC’s on the first stage. Unlike 
the two-dimensional caserthe standard treatment of 
the external boundary in three dimensions is based 
on merely the locally one-dimensional characteristics 
analysis and extrapolation (the point-vortex model 
is not applicable). 

In Table IT we present the numerical results ob- 
tained with the two types of ABC’s for the subcrit- 
ical flow around the ONERA M6 wing. As one can 
clearly seeMie DPM-based ABC’s in three space di- 
mensions outperform the standard approach as they 
do in two dimensions: the new boundary conditions 
are capable of producing accurate results on very 
small computational domainsr whereas the perfor- 
mance of the standard technique deteriorates as the 
domain shrinks. Analogously to the two-dimensional 
caserthe three-dimensional DPM-based ABC’s can 
be calculated for any shape of T by means of the 
same procedure (universality) and can be easily com- 
bined with the existing solver (TLIS3D). Let us ad- 
ditionally note that for this series of computationsr 
the dimensions of all three grids are the samergrids 
of smaller extent are obtained by scaling down the 
original 10 chords grid. The concentration of nodes 
near the wing that results from the down-scaling ob- 
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viously contributes to the general improvement in 
the accuracy of the solution. 

We have also conducted some three-dimensional 
computations for the standard transonic case for 
ONERA M6 wing: M 0 = 0.84; Re 0 = 11.7 • 10 6 ; 
a = 3.06°. Here we could not bring the artificial 
boundary as close to the solid surface as in the fore- 
going subsonic case (see Table 1) because our far- 
held treatment is purely subsonic and the boundary 
cannot intersect the supercritical zone. We there- 
fore calculated the solution for two computational 
domains of the average radii of approximately 10 
and 3 root chords of the wingr respectively. The 
results summarized in Table 2 clearly demonstrate 
that for the small computational domains the DPM- 
based ABC’s generate much more accurate solutions 
that the standard boundary conditions do. Noterin 
this case the grids for the different domains have 
different dimensionsrand the smaller 197 x 49 x 33 
grid (3 chords) is now an exact subset of the big- 
ger 203 x 57 x 33 grid (10 chords). This is done in 
order to eliminate any influence that the change of 
the grid in the near field could possibly exert on the 
solution. 

Besides the improvement of accuracyr the appli- 
cation of the DPM-based ABC’s to transonic flow 
computations on the small (3 chords) computa- 
tional domain yielded much higher convergence rate 
of the residual (continuity equation )T as well as 
much faster convergence of other quant. itiesrinclud- 
ing those deemed as sensitivePe.g.rthe number of su- 
personic points in the domain. In Figures 3a. and 3br 
we show the convergence history for this supercrit- 
ical flow variant.. One can see that the convergence 
for the standard boundary conditions is poor; there- 
fore the corresponding force coefficients in Table 2 
are given with the error bands indicated. 

For the 10 chords doma.inrt.he DPM-ba.sed ABC’s 
also provide for some convergence speedupPalthough 
the difference between the two ABC’s techniques is 
less dramatic here. This is reasonable because one 
could generally expect, that, the bigger the computa- 
tional doma.inrt.he smaller is the influence that, the 
external boundary conditions exert, on the numeri- 
cal procedure. The convergence history for the 10 
chords computations is shown in Figures 4a. and 4b. 

Not.erfrom Figure 3b one can conclude that, on the 
small domain the: two algorithms converge to quite 
different, solut.ionsr whereas Figure 4b allows one to 
assumeThat. on the big domain the final solutions are 
close to one another. The data, from Table 2 corrob- 
orate these conclusions. This behavior again fits into 
the aforementioned concept, that, the impact, of the 
ABC’s decreases as the domain size increases. In the 


fut.urer we plan on running some transonic compu- 
tations for the domain bigger than 10 root, chords of 
the wing. This may provide for an experimental jus- 
tification of the statement, that, the solutions for the 
two different, types of ABC’s approach one another 
as the domain enlarges. 



Work 

Figure 3a. ONERA M6: M 0 = 0MTRe o = llJ-lO'T 
a = 3.06°. Convergence history for the residual of 
the continuity equation. Average domain “radius” 
is 3 root, chords of the wing; grid 197 x 49 x 33. 



Work 

Figure 3b. ONERA M6: M 0 = 0MTRe o = 11.7 • 
lO'To- = 3.06°. Convergence history for the number 
of supersonic nodes in the domain. Average domain 
“radius” is 3 root, chords of the wing; grid 197 x.49 x 
33. 

Miscellaneous Issues. Herer we will primarily 
discuss the computational cost, of the DPM-ba.sed 
ABC’s and possible ways for its reduction. 

The cost, of the DPM-ba.sed ABC’s in two space 
dimensions is modest.; the ABC’s may add about. 
10%) to the cost, of the original integration proce- 
dure. In three space dimensionsrt.he relative cost, of 
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Figure 4a. ONERA M6: M 0 = 0MTRe o = 11.7-10'T 
a = 3.06°. Convergence history for the residual of 
the continuity equation. Averag# domain “radius” 
is 10 root chords of the wing; grid 209 x 57 x 33. 
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Work 

Figure 4b. ONERA M6: M 0 = 0MTRe o = 11.7 • 
lO'To- = 3.06°. Convergence history for the number 
of supersonic nodes in the domain. Average domain 
“radius” is 10 root chords of the wing; grid 209 x 
57 x 33. 


the DPM-based ABC’s is so far higher; they typi- 
cally add from 25%) to 30%) of extra CPU time. In 
many casesrhoweveiTthis increase can be compen- 
sated for by the convergence a.ccelera.t.ionr see Fig- 
ures 3 and 4. MoreoveiTwe expect that some changes 
to the ABC’s algorithm can still be done so that the 
vect.oriza.t.ion is exploited more effectively and there- 
fore the boundary operators are computed faster. 
Apart from the algorithm changes there are some 
“indirect” ways for reducing the computational cost, 
of the DPM-based ABC’s. 

Parallel implementation of the ABC’s on the 
multi-processor machines is a relatively easy task 


because the solution algorithm for the linear AP is 
readily parallelizable. For two space dimensionsFan 
up to five times speedup compared to the single- 
processor version has been achieved on an eight- 
processor CRAY Y-MP. 

In three space dimensionsr carrying out multiple 
computations on the same grid is very often the case 
(many different flow regimes in one geometric set- 
ting). In so doingrt.he external geometry that in- 
fluences the ABC’s is fixed and the boundary con- 
ditions depend only on the aerodynamic para.me- 
t.ersr e.gT Mach number Mo and the angle of at- 
tack a. Thereforer the entry- wise interpolation of 
the matrices of boundary operators T with respect 
to these parameters has a substantial promise from 
the standpoint of cost, reduction. Indeedr after a. 
noticeable startup expense for calculating the “ref- 
erence” T’s for some selected values of the pa.ra.me- 
t.ersFa.ny other matrix needed for any specific regime 
within the prescribed range of the parameters can 
be obtained for virtually no extra, cost, by means of 
the interpolation. This methodology has been tested 
for the two-dimensional supercritical flow around the 
RAE2822 airfoil . The results seem very promis- 
ing: for the small computational domain the accu- 
racy provided by the operators T interpolated with 
respect, t.o Mo and a is only slightly worse than the 
accuracy obtained using genuine DPM-based ABC’s 
(i.e.rt.he operator calculated for the given specific 
values of the parameters) and is still much better 
than the accuracy obtained on the basis of the point- 
vort.ex model. 

Fina.llyrwe mention that, the nonloca.lit.y of the 
DPM-based ABC’s is expressed by the fact, that the 
matrix T (see (10)) is dense. Howe veiT this ma- 
trix is typically structured. From the standpoint, 
of physicsr the structure of T reflects the, simple 
consideration that, each specific node influences its 
close neighbors stronger than it. influences the re- 
mote points. For system (4a.)rt.he operator T from 
(10) is composed of 5 x 5 blocks as shown in Fig- 
ure 5 (provided that, the vectors of variables along 
the boundaries T and Ti are arranged properly: first, 
p’s for all the nodesrt.hen w’sret.c.ra.ndrHna.llyrp’s). 
The qualitative meaning of each block in Figure 5 is 
how one variable of the corresponding pair influences 
another one along the entire boundary. Stronger in- 
terdependence between closer nodes results in the 
fact. that, the near-diagonal entries for each block 
are considerably larger than the off-diagonal ones (ifT 
a.dditiona.llyrthe order of nodes for T and T i is the 
same). Thb latter observation ha.srin factTbeen cor- 
roborated computationally (for the two-dimensional 
case). 


9 


American Institute of Aeronautics and Astronautics 



p~p 

u ~ P 

v ~ P 

w ~ p 

p~ p 

p~u 

u ~ u 

V ~ u 

w ~ u 

p~u 

p~v 

U ~ V 

V ~ V 

W ~ V 

p~v 

p ~ w 

u ~ w 

V ~ w 

w ~ w 

w ~ p 

p~ p 

u~p 

v~p 

w ~ p 

p~p 


Figure 5. Qualitative block structure of the bound- 
ary operator T for three-dimensional flow computa- 
tions. 

To effectively calculate matrix- vector products for 
the structured matrices of the foregoing typer one 
can use multiresolution-based data compression a.l- 
gorithmsr e.g.T the one proposed by Harten and 

OA 

Yad-Shalom . In Ref. 19r we applied this algo- 
rithm to the model example of the matrix T for 
two-dimensional scalar Yukawa equation. As could 
be seen from the nested multiresolution representa- 
tion for the model operator T (see Ref. 19)Fa.ll 
the nonzero entries are concentrated near the main 
diagonal of each nested block; as we recede the diag- 
onal the corresponding values rapidly fall below the 
double-precision machine accuracy (10 -13 ) thresh- 
old. Therefore rt.his nested structure can be used 
for the efficient calculation of matrix- vector products 
with consecutively improving resolution. The con- 
centration of nonzero entries near the main diagonal 
also implies effective localization of the boundary 
conditions. This feature may have certain promise 
for the future use in the view of the multi-block grids. 

The next major challenge from the standpoint 
of treating numerically the infinite-domain formula- 
tions is presented by the time-dependent, problems. 
As mentioned in the beginningrt.he issue of ABC’s 
for the time-dependent, problems is significant, not. 
only in CFD but. also in other areas of scientific com- 
put.ingr e.gT those originating from acoustics and 
electrodynamics. 

The primary difficulty that, distinguishes the time- 
dependent. problems from the st.ea.dy-st.a.t.e ones is 
that, here the exact. ABC’s are nonlocal not. only in 
space but. also in time (except, in the simplest, one- 
dimensional problems). As concerns the approxi- 
mate a.pproa.chesrt.he same basic trend that, is rele- 
vant. to the st.ea.dy-st.a.t.e problems also holds for t.hft 


time-dependent, ones: more accurate boundary pro- 
cedure requires more of the nonlocal nature of the 
ABC’s to be taken into account.. 

Of courser the nonloca.lit.y of the ABC’s in time 
presents a. very serious computational obst.a.clerbe- 
ca.use as the solution evolves such boundary condi- 
tions would become more and more expensive from 
the standpoints of both memory and computer time 
requirements. Thereforer the most, crucial numeri- 
cal issue in constructing the time-dependent. ABC’s 
is how to effectively restrict, the nonloca.lit.y of the 
boundary conditions in time. This is an interesting 
subject, for the future research. 

4 Conclusions. 

The new global ABC’s for calculating st.ea.dy-st.a.t.e 
external viscous flows in three space dimensions have 
been constructed on the basis of the difference poten- 
tials method. The approach generalizes and extends 
our previous two-dimensional results. 

The new ABC’s are capable of greatly reducing 
the size of the computational domain (compared to 
the standard methods) while still maintaining high 
accuracy of the numerical solution. This size re- 
duction amounts to either the possibility of refining 
the grid in the near fieldr which potentially leads 
to the improvement, in accuracy or usage of the 
smaller- dimension grids without, compromising the 
accuracy. MoreoveiTt.he DPM-ba.sed ABC’s may no- 
ticeably improve the convergence rate of the multi- 
grid iteration procedure. Fina.llyrt.he new boundary 
conditions appear geometrically universal and easy 
to incorporate in the structure of the existing flow 
solvers. The properties of the new ABC’s have been 
corroborated experimentally by computing the sub- 
sonic and transonic flows past, the ONERA M6 wing 
using the NASA-developed code TLIS3D. 

A more detailed description of the theoretical 
foundations of the DPM-ba.sed a.lgorit.hmra.s well as 
a. wider selection of the computational result.sr will 
be presented in the forthcoming paper Ref. 25. 
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